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ABSTRACT 

We report on a numerical study of viscous fluid accretion onto a black hole. The flow is axisymmetric 
and uses a pseudo-Newtonian potential to model relativistic effects near the event horizon. The numerical 
method is a variant of the ZEUS code. As a test of our numerical scheme, we are able to reproduce results 
from earlier, similar work by Igumenshchev and Abramowicz and Stone et al. We consider models in 
which mass is injected onto the grid as well as models in which an initial equilibrium torus is accreted. In 
each model we measure three "eigenvalues" of the flow: the accretion rate of mass, angular momentum, 
and energy. We find that the eigenvalues are sensitive to r£ n , the location of the inner radial boundary. 
Only when the flow is always supersonic on the inner boundary are the eigenvalues insensitive to small 
changes in rj n . We also report on the sensitivity of the results to other numerical parameters. 

Subject headings: accretion disks, black hole physics, hydrodynamics, turbulence, galaxies: active 



1. INTRODUCTION 

Black hole accretion flows are the most likely cen- 
tral engine for quasars and active galactic nuclei (AGN) 
(Zel'dovich 1964; Salpeter 1964). As such they are the sub- 
ject of intense astrophysical interest and speculation. Re- 
cent observations from XMM-Newton, Chandra, Hubble, 
VLBA, and other ground- and space-based observatories 
have expanded our understanding of the time variability, 
spectra, and spatial structure of AGN. Radio interferome- 
try, in particular, has been able to probe within a few hun- 
dred gravitational radii (GM/c 2 ) of the central black hole, 
e.g. Lo et al. (1998); Junor et al. (1999); Doeleman et al. 
(2001). Despite these observational advances, only instru- 
ments now in the concept phase will have sufficient angu- 
lar resolution to spatially resolve the inner accretion disk 
(Rees 2001). And so there remain fundamental questions 
that we can only answer by folding observations through 
models of AGN structure. 

All black hole accretion flow models require that angu- 
lar momentum be removed from the flow in some way so 
that material can flow inwards. In one group of models, 
angular momentum is removed directly from the inflow 
by, e.g., a magneto-centrifugal wind (Blandford & Payne 
1982). Here we will focus on the other group of models in 
which angular momentum is diffused outward through the 
accretion flow. 

It has long been suspected that the diffusion of angular 
momentum through an accretion flow is driven by turbu- 
lence. The a model (Shakura & Sunyaev 1973) introduced 
a phenomenological shear stress into the equations of mo- 
tion to model the effects of this turbulence. This shear 
stress is proportional to aP, where a is a dimensionless 
constant and P is the (gas or gas + radiation) pressure. 
This shear stress permits an exchange of angular momen- 
tum between neighboring, differentially rotating layers in 
an accretion disk. In this sense it is analogous to a viscos- 
ity (see also Lynden-Bell & Pringle (1974)) and is often 
referred to as the "anomalous viscosity." 
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The a model artfully avoids the question of the origin 
and nature of turbulence in accretion disks. This allows 
useful estimates to be made absent the solution to a diffi- 
cult, perhaps intractable, problem. Recently, however, sig- 
nificant progress has been made in understanding the ori- 
gin of turbulence in accretion flows. It is now known that, 
in the magnetohydrodynamic (MHD) approximation, an 
accreting, differentially rotating plasma is destabilized by 
a weak magnetic field (Balbus & Hawley 1991; Hawley & 
Balbus 1991). This magneto-rotational instability (MRI) 
generates angular momentum transport under a broad 
range of conditions. Numerical work has shown that in 
a plasma that is fully ionized, which is likely the case for 
the inner regions of most black hole accretion flows, the 
MRI is capable of sustaining turbulence in the nonlinear 
regime (Hawley & Balbus 1991; Hawley et al. 1995; Haw- 
ley 2000; Hawley & Krolik 2001). 

Studies of unmagnetized disks have greatly reduced the 
probability that a linear or nonlinear hydrodynamic in- 
stability drives disk turbulence. While there are known 
global hydrodynamic instabilities that could in principle 
initiate turbulence, these have turned out to saturate at 
low levels or require conditions that are not relevant to an 
accretion disk near a black hole. As of this writing, no 
local, linear or nonlinear hydrodynamic instabilities that 
transport angular momentum outwards are known to exist 
in Keplerian disks (Balbus & Hawley 1998). 

Work on magnetized disks has now turned to global nu- 
merical models. These are possible thanks to advances in 
computer hardware and algorithms. Recent work by Haw- 
ley (2000, 2001), Stone & Pringle (2001), and Hawley & 
Krolik (2001) considers the evolution of inviscid, nonrela- 
tivistic MHD accretion flows in two or three dimensions. 
Some of this work uses a pseudo-Newtonian, or Paczyn- 
ski & Wiita (1980), potential as a model for the effects of 
strong-field gravity near the event horizon. 

Other work on global models has considered the equa- 
tions of viscous, compressible fluid dynamics as a model for 
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the accreting plasma (Igumenshchev & Abramowicz 1999; 
Stone et al. 1999; Igumenshchev & Abramowicz 2000; Igu- 
menshchev et al. 2000). The viscosity is meant to mock 
up the effect of small scale turbulence, presumably gener- 
ated by magnetic fields, on the large scale flow. In light 
of work on numerical MHD models, this may seem like 
a step backwards. The MHD models, however, are com- 
putationally expensive and introduce new problems with 
respect to initial and boundary conditions. It therefore 
seems reasonable to investigate the less expensive a based 
viscosity models. In this paper we investigate axisymmet- 
ric, numerical, viscous inflow models. 

This work was motivated by the earlier work of Igu- 
menshchev & Abramowicz (1999, 2000) and Stone et al. 

(1999) , hereafter referred to as IA99, IA00, and SPB99, 
respectively (IA99 and IA00 are collectively referred to as 
IA, in which case SPB99 is simply referred to as SPB). 
These authors studied similar viscous inflow models yet 
found different radial scaling laws for radial velocity, den- 
sity, and angular momentum. They also found different 
values for the accretion rate of mass and angular momen- 
tum. They used different experimental designs, however. 
We set out to discover whether the results from these au- 
thors differed due to numerical methods or model param- 
eters. 

Along the way, we took a systematic approach to study- 
ing numerical parameters and boundary conditions. One 
particular point of concern, which will be described in 
greater detail below, is the inner boundary condition. This 
lies in the energetically dominant portion of the flow, so er- 
rors there can potentially corrupt the entire model. In this 
paper we show that aspects of results presented by other 
researchers are sensitive to model and numerical parame- 
ters. These results should be useful to others contemplat- 
ing large-scale numerical models of black hole accretion. 

The paper is organized as follows. In §2 we discuss our 
models. In §3 we discuss numerical methods. In §4 we dis- 
cuss a fiducial solution and results from a survey of other 
solutions. In §5 we summarize our results. 

2. MODEL 

We are interested in modeling the plasma within a few 
hundred GM/c 2 of a black hole. We will consider only 
axisymmetric models (the work of Igumenshchev et al. 

(2000) suggests that 2D and 3D viscous models give sim- 
ilar results). Throughout we use standard spherical polar 
coordinates r, 9, and <p. 

We solve numerically the axisymmetric, nonrelativistic 
equations of compressible hydrodynamics in the presence 
of an anomalous stress II, which is meant to model the 
effects of small-scale turbulence on the mean flow. The 
governing equations then express the conservation of mass 
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momentum, 



and energy, 



P ^ = -vp - pv* - v • n, 
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-(P + u)(V-v) 



(1) 
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(3) 



Here, as usual, D/Dt = d/dt + v • V is the Lagrangian 
time derivative, p is the mass density, v is the velocity, 



u is the internal energy density, P is the pressure, and 
\P is the gravitational potential. The dissipation function 
$ is given by the product of the anomalous stress tensor 
II with the rate-of-strain tensor e (given explicitly in the 
Appendix) 



(4) 



(sum over indices) where the anomalous stress tensor is 
the term-by-term product 



-2pveijS, 



(5) 

(no sum over indices) where SV, is a symmetric matrix 
filled with or 1 that serves as a switch for each compo- 
nent of the anomalous stress. The equation of state is 



P= (7-l)u. 



(6) 



For the gravitational potential we use the pseudo- 
Newtonian potential of Paczynski & Wiita (1980): = 
-GM/(r-r g ) (here r g = 2GM/c 2 ). This potential repro- 
duces features of the orbital structure of a Schwarzschild 
spacetime, including an innermost stable circular orbit 
(ISCO) located at r = 6GM/c 2 . In a few cases we use the 
Newtonian potential W = —GM/r for comparison with 
others' work (IA and SPB exclusively use a Newtonian 
potential). 

We must now make some choices for the anomalous 
stress tensor. One might argue on very general grounds 
for a Navier-Stokes prescription, and indeed IA and we 
use a prescription where all elements of S are 1 (the "IA 
prescription"). SPB, on the other hand, use the Navier- 
Stokes prescription with all components zero except S r ^, 
Sec/,, S^r, and S^g (the "SPB prescription"). SPB justify 
this choice by arguing that it more appropriately models 
MHD turbulence; this was later supported by results pre- 
sented in Stone & Pringlc (2001). 

We must also choose a viscosity coefficient. We consider 
three different prescriptions: one similar to those chosen 
by IA; a second viscosity coefficient similar to that chosen 
by SPB; and a third, similar form that vanishes rapidly 
near the poles. Explicitly, 



v = a(c 2 /Q, K ), 
v = a(p/p )n a r 2 , 

v = a(c 2 /n K )sm 3 / 2 (6), 



(7) 
(8) 
(9) 



are the IA, SPB, and MG prescriptions, respectively, 
where c s = \/-fP/p is the sound speed, and 



O 2 = 
n K — 



GM 



i a* 



r dr r(r — r g ) 2 



(10) 



is the "Keplerian" angular velocity. Here po and fio are 
values of p and f2 at a fiducial radius r . 

The choice of viscosity coefficient for IA and MG is 
based, as usual, on dimensional arguments. SPB99's 
choice focuses the viscosity where most of the matter is, 
a numerical convenience. Our MG prescription is a small 
modification of the IA prescription to concentrate the vis- 
cosity toward the equator. These choices are to a large 
extent arbitrary, although one might attempt to motivate 
the choice by comparison with MHD simulations, as do 
Stone & Pringle (2001). Nevertheless, some dynamical 
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properties of MHD turbulence, such as the elastic prop- 
erties that produce magnetic tension and hence Alfven 
waves, can never be modeled with a viscosity. 

The model has boundaries at 6 = ±n/2 and at r = 
Tin ! r out ■ At the 9 boundaries we use the usual polar 
axis boundary conditions. At the radial boundaries we 
use "outflow" boundary conditions; ideally these boundary 
conditions should be completely transparent to outgoing 
waves. 

Fuel for the accretion flow must be provided cither in 
the initial conditions or continuously over the model evo- 
lution. Some global numerical accretion flow models have 
started with an equilibrium torus. Examples include Haw- 
ley (1991), Hawley et al. (1995), Hawley (2000), Haw- 
ley (2001), Hawley et al. (2001), and Hawley & Krolik 
(2001). Others have started with an initial configuration 
of matter that is not in equilibrium. For example, mat- 
ter may be placed in orbit about the black hole, but with 
sub-Keplerian angular momentum, so that once the sim- 
ulation commences it immediately falls toward the hole. 
Examples of this approach include Hawley et al. (1984), 
Koidc et al. (1999), Koide et al. (2000), and Meier et al. 
(2001). This approach may enhance transients associated 
with the choice of initial conditions, although it can also 
be physically well motivated, as in studies of core-collapse 
supernovae. One can also inject fluid continuously onto 
the computational grid over the course of the evolution. 
Examples of this include IA99 and IA00. One might also 
use an inflow boundary condition at the outer radial edge 
of the grid, as in Blondin et al. (2001). The main advan- 
tage of injection models is that they allow one to achieve a 
steady, or statistically steady, state. In this paper we will 
consider only the equilibrium tori and on-grid injection 
models. 

The equilibrium tori (Papaloizou & Pringle 1984) are 
steady-state solutions to the equations of inviscid hydro- 
dynamics. They assume a polytropic distribution of mass 
and internal energy, P = Kp 1 , and a power-law rota- 
tion profile, cx (rsin#)~ 9 . The radial and meridional 
components of the velocity vanish. There are 5 parame- 
ters that describe the torus: (1) the location of the torus 
pressure maximum ro ; (2) the location of the innermost 
edge of the torus, r tjin < r ; (3) the maximum value of 
the density, p — p(r ); (4) the angular velocity gradient 
q = d\nQ/d\nR; and (5) the entropy constant K. 

On-grid injection adds matter to the model at a con- 
stant rate. The matter is injected with a non-zero spe- 
cific angular momentum and with zero radial or meridional 
momentum in a steady pattern p(r, 9) which is typically 
symmetric about the equator. Parameters for this scheme 
include: (1) a characteristic radius for injection r in j] (2) 
the rate of mass injection M in j; (3) the specific angular 
momentum of the injected fluid, v$ = JitVLk- We usu- 
ally set fi = 0.95, so that the fluid circularizes near Ti n j. 
This restricts transients associated with circularization to 
the outer portions of the computational domain; (4) the 
internal energy of the injected fluid, u = f^p^ . We always 
set fi = 0.2 so that the fluid is marginally bound, i.e. has 
Bernoulli parameter Be = (l/2)v 2 + c 2 s /(-y — 1) + ^ < 0. 

One must also choose the injection pattern p(r, 9). IA99 
choose a radially narrow region, but do not explicitly give 
p(r, 9). We use a Gaussian, but found that none of the re- 



sults are sensitive to the precise profile. The accretion rate 
of mass, energy, and angular momentum are insensitive to 
large changes in the size of the injection region except for 
the extreme cases of filling the entire width or injecting 
in 2 locations. These extreme cases arc sufficiently differ- 
ent to be referred to as a completely different model; they 
lead to a qualitative change in the flow. For example, a 
full range 9 injection region has matter that will collide 
with any outflow at the poles. A bipolar injection leads 
to an equatorial outflow. Our models have radial width 
oy = 0.05(r in + r out )/2 and a = w/8. 

3. NUMERICAL METHODS 

Our numerical method is based on ZEUS-2D (Stone & 
Norman 1992) with the addition of an explicit scheme for 
the viscosity. ZEUS is an operator-split, finite-difference 
algorithm on a staggered mesh that uses an artificial vis- 
cosity to capture shocks (in addition to the anomalous 
viscosity in equations [2] ) . This algorithm guarantees that 
momentum and mass are conserved to machine precision. 
Total energy is conserved only to truncation error, so total 
energy conservation is useful in assessing the accuracy of 
the evolution. 

The inner and outer radial boundary conditions are im- 
plemented by copying primitive variable values (p, u, and 
v) from the last zone on the grid into a set of "ghost 
zones" immediately outside the grid. Inflow from outside 
the grid is forbidden; we set v r (ri n ) — if v r (ri n ) > and 
v r (r ut) = if v r (r out ) < 0. Since we expect inflow on the 
inner boundary, this switch should seldom be activated. 
We have found that frequent activation of the switch is 
usually an indication of a numerical problem. 

We use a radial grid uniform in log(r — r g ). We require 
that dr{r)/(r in — r g ) < 1/4 so that the structure of the 
pseudo-Newtonian potential is well resolved. The grid is 
uniform in 9. The grid has N r x Ng zones. 

3.1. Numerical Treatment of Low Density Regions 

Like many schemes for numerical hydrodynamics, ZEUS 
can tolerate only a limited dynamic range in density. It 
is therefore necessary to impose a density minimum pfi 
to avoid small or negative densities. Our procedure for 
imposing the floor is equivalent to adding a small amount 
of mass to the grid every time the floor is invoked. Mass 
is added in such a way that momentum is conserved. To 
monitor the effect of the density floor, we track the rate of 
change of total mass and total energy (from kinetic energy 
change) due to this procedure, Mfi and Efi. 

We set p f i = 10" 10 M inj c 3 /(GM) 2 for injection runs 
and pfi = f0~ 5 /?o for torus runs. Lower values for pfi do 
not lead to a significant change in the solution. Larger 
values of pfi give Mfi ~ M, the accretion rate through 
the inner boundary. The atmosphere also becomes more 
massive and begins to affect torus stability- vertical oscil- 
lations are excited in the inner disk by a Kelvin-Helmholtz 
like instability. This should be avoided. 

We must also surround the torus in a low density atmo- 
sphere in the initial conditions. The density of the atmo- 
sphere is pfi and the internal energy density is u ~ U p^ , 
where U is a constant fraction of order unity (e.g. IA and 
we choose U — 0.2). The addition of the atmosphere has 
no effect on the solution since the mass source's evolution 
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eventually dominates the flow everywhere. SPB99 choose 
a different method of constructing the initial atmosphere 
but obtain late-time results that are similar to ours. 

ft is also necessary to impose a floor u fi on the internal 
energy density. This we take to be the minimum value of 
u in the initial atmosphere. As for the mass, we track the 
rate of change of total energy due to the internal energy 
floor, that along with the kinetic energy is included in Efi. 

3.2. Diagnostics 

Global numerical simulations of accretion flows are com- 
plicated; it is possible to measure many quantities associ- 
ated with the flow. Some are astrophysically relevant, and 
some are not. In our view particular interest attaches to 
the time-averaged flux of mass, energy, and angular mo- 
mentum through the inner boundary. Physically, these are 
directly related to the luminosity of the accretion flow and 
the rate of change of mass and angular momentum of the 
central black hole. As described by Narayan & Popham 
(1993), these are in a sense the nonlinear "eigenvalues" of 
the model. 

The mass accretion rate is 

M = J pv-dS, (11) 

s 

where S is the inner radial surface of the computational 
domain. The total energy accretion rate is 

E = J ((^v 2 +h+ <F)pv + n • v) • dS, (12) 

s 

where h = (u + P)/p = ju/p is the specific enthalpy with 
our equation of state. The angular momentum accretion 
rate is 

L = Jrsme(pvv (t> +U-4>)-dS. (13) 

s 

It is also sometimes useful to focus on the reduced eigen- 
values I = L/M and e = E/M. These value of mass, 
energy, and angular momentum are recorded at about 2 
grid zones away from R in . This avoids any error that may 
occur when evaluating directly on the boundary where the 
inflow boundary condition is applied. 

We also track volume-integrated quantities, the flux of 
mass, energy, and angular momentum across all bound- 
aries, and floor added quantities in order to evaluate the 
consistency of the results. Mass and angular momentum 
are conserved to machine precision, although "machine 
precision" implies a surprisingly large random walk in the 
integrated quantities over the full integration because the 
calculation requires millions of timesteps. 

Total energy is conserved to truncation error, not ma- 
chine precision, and thus is a useful check on the quality 
of the simulation. Total energy conservation implies 

E err — E vo i + E + E out — Efi, (14) 

where E vo i is the rate of change of the volume integrated 
total energy, E out is the flux of total energy through the 
outer radial boundary, and Efi is the rate of total energy 
added due to the kinetic energy change (because of the 
mass density floor) and internal energy density floor. Ide- 
ally, E err — 0. Truncation errors can (and do) lead to 



cumulative, rather than random, changes in the total en- 
ergy. A useful gauge of the magnitude of these errors is 
E err /E. For all runs we performed the error rate is within 
10% of W~ 5 c 3 /GM for a torus run and within 10% of 
10~ 4 c 3 /GM for a viscous injection run. 

3.3. Code Tests 

Our version of ZEUS reproduces all hydrodynamic test 
results from Stone & Norman (1992), including their 
spherical advection and Sod shock tests. We also find 
excellent agreement with steady spherical accretion solu- 
tions, i.e. Bondi flow /citcpbondi52. An inviscid equi- 
librium torus run also persists for many dynamical times 
with insignificant deviations from the initial conditions. 

We have parallelized our code using the MPI message 
passing library. On the Origin 2000 at NCSA we are able 
to achieve about 2.5 x 10 zone updates per second using 
240 CPUs, or about 35 GFLOPs. This is 159 times faster 
than the single CPU speed, which represents a parallel 
efficiency of 66%. 

4. RESULTS 

The initial motivation for undertaking this calculation 
was to understand differences between results reported in 
IA and SPB. Using our code, which is based on the same 
algorithm used in SPB's calculations, we ran a series of 
tests attempting to reproduce SPB's results. These test 
calculations used all of SPB's model choices, including 
SPB's viscosity prescription, a Newtonian potential, and 
a torus for the mass source. We were able to reproduce 
most quantitative results reported in SPB99's torus cal- 
culations. This includes their radial scaling laws. For ex- 
ample, in a model that is identical to SPB99's Run B, we 
find Mar, p oc r°, c 2 s oc r -1 , oc r -1 / 2 , and \v r \ oc r _1 . 
These power law slopes are identical to those reported by 
SPB99. As another example, we found M = 1.23 x 10~ 3 
torus masses per torus orbit at the pressure maximum for 
a model identical to SPB's Model A (their fiducial model); 
SPB report M = 1.0 x 10~ 3 in the same units. Given the 
fluctuations in mass accretion rate, our value and SPB's 
value are fully consistent. We were even able to reproduce 
certain numerical artifacts associated with the inner radial 
boundary, such as a density drop and temperature spike 
near the inner boundary. 

Recall that SPB evolve an initial torus and allow it to 
accrete; IA use a different experimental design in which 
matter is steadily injected onto the grid. They also use 
a different viscosity prescription. We ran a second series 
of test calculations attempting to reproduce IA's results. 
These test calculations used all of IA's model choices, in- 
cluding viscosity prescription, Newtonian potential, etc. 
We were able to reproduce all of IA99's calculations ex- 
cept those that include thermal conduction (which we did 
not attempt to reproduce). In each case we found that the 
qualitative nature of the flow is similar to that described in 
IA99. In particular, we agree on which models are stable 
and unstable and which models exhibit outflows. We also 
find qualitative agreement with their contour plots of, e.g., 
density pressure, mass flux, and Mach number. We also 
find qualitative agreement with their radial run of c s /Vk 
and specific angular momentum. Our results do not agree 
precisely, but this is likely due to small differences in mass 
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injection scheme (because IA99 do not give their p(r,8)). 
Finally, we can also reproduce the radial scalings given in 
IA00 for their model A. 

The most significant difference between the results of IA 
and SPB was due to the choice of anomalous stress pre- 
scription, as might have been anticipated. Qualitatively, 
the stress components that are included in IA and not 
SPB tend to smooth the flow and suppress turbulence. 
Thus SPB99's simulations result in more vigorous convec- 
tion than simulations performed by IA. The choice of mass 
supply (torus vs. injection) also leads to a significant dif- 
ference between IA and SPB's results; this is discussed in 
more detail below. 

The fact that we can reproduce both SPB's and IA's re- 
sults using a single code is consistent with the hypothesis 
that differences between their reported results (e.g. the 
lower degree of convection reported in IA than SPB, and 
the differences in radial scaling laws) is due to differences in 
experimental design and viscosity prescription rather than 
numerical methods. While we cannot completely rule out 
the possibility that SPB, IA, and we have made identical 
experimental errors, this seems unlikely. This compari- 
son thus lends credibility to SPB, IA, and our numerical 
results. 

4.1. Fiducial Model Evolution 

We now turn from reproducing earlier viscosity mod- 
els to considering new aspects of our own models. First, 
consider the evolution of a "fiducial" model (Run A in 
Table Al and Table A2). The fiducial model has r in = 
2.7GM/c 2 , r out = 600GM/c 2 , r inj = 495GM/c 2 , 7 = 3/2, 
a = 0.1, N r — 108, Ng = 50. It uses a pseudo- Newtonian 
potential, mass is supplied by injection, and the viscos- 
ity prescription follows IA. It was run from t = to 
t '= 7.3 x 10 5 GM/c 3 . 

Run A is similar to IA99's "Model 5" , except that it uses 
a pseudo-Newtonian potential. In a statistically steady 
state the flow is characterized by a quasi-periodic out- 
flow. Hot bubbles form at the interface between bound 
(Bernoulli parameter Be = (l/2)v 2 + c 2 J(-f - 1) + * < 0) 
and unbound (Be > 0) material. These hot bubbles are 
buoyant and move away from the black hole. This appears 
to be a low-frequency, low wavenumber convective mode 
(IA refer to it as a "unipolar outflow" ) . Higher wavenum- 
ber convective modes are evidently suppressed by the vis- 
cosity. 

Figure Al shows time-averaged plots of various quanti- 
ties in the fiducial run. The time average is performed from 
200 equally spaced data dumps from t = 4.3 x 10 5 GM/c 3 
to t = 7.3 x 10 5 GM/c 3 . We show only the region r in < 
r < 30GM/c 2 , whereas the computational domain is much 
larger: r„ < r < 600GM/c 2 . The injection point is lo- 
cated far outside the plotted domain at r = 496GM/c 2 . 
Because of the strong time-dependence of the flow in the 
fiducial run, the flow at any instant may look very different 
from these time averaged plots. 

The upper left panel in Figure 1 shows the average den- 
sity; notice that the density is not symmetric about the 
equator. This is because the flow involves long-timescale 
quasi-periodic variations which arc not quite averaged out 
over the course of the run. The upper right corner shows 

3 Averaging over \9 — tt\ < 7r/36 produces nearly identical results, but 



the Bernoulli parameter Be. Dotted lines are negative; 
notice that there is a substantial amount of fluid near the 
equator that is unbound in the sense that Be > 0. Never- 
theless, this material is still flowing inward in a nearly lam- 
inar fashion. Near the poles, the time-averaged Be < 0, 
but this region experiences large fluctuations. Polar out- 
flows arc typically associated with positive fluctuations 
in Be. The lower left panel shows the scaled mass flux 
r 2 s'm8(pv). Notice that much of the mass flux is along 
the surfaces of the inflow rather than along the equator. 
The lower right panel shows the scaled angular momentum 
flux r 3 sin 2 9(pwv t j 1 + II • </>). As for the mass flux, most of 
the activity is along the surface of the flow. 

Figure A2 shows the time series of the reduced eigen- 
values: M/Minj, e = E/(Mc 2 ), and I = Lc/(GMM). 
Also shown as dashed lines are the thin disk values for 
e and I. These assume a thin, cold disk terminating at 
r = 6GM/c 2 . The low value of I is due to two effects. 
First, the disk is already sub-Keplerian by the time the 
flow reaches the innermost stable circular orbit. In ad- 
dition, there are residual viscous torques in the plunging 
region that lower the specific angular momentum of the 
accreted material (see Figure A3, below). Notice that Fig- 
ure A2 shows a smooth evolution that varies on a timcscale 
t « 4 x 10 4 at late time. The largest variations in mass 
accretion rate are related to the appearance of large con- 
vective bubbles. 

Figure A3 shows the 8 and time averaged run of sev- 
eral quantities with radius. The averages are taken over 
4.3xl0 5 GM/c 3 < t < 7.3xl0 5 GM/c 3 and \8-w/2\ < ir/6 
. The upper left plot shows the run of density. Notice 
that here, as for the other quantities, there is a spike near 
ri n j, an intermediate region, and then an inner, roughly 
power-law region. The upper right plot shows (c s /c) 2 ; the 
lower left shows |w r |/c. Notice that the radial velocity ex- 
ceeds the speed of light at the inner boundary. Similarly 
the azimuthal velocity v^/c shown in the lower right panel 
approaches the speed of light. Also shown in that panel is 
the circular velocity (dashed line). Evidently the flow is 
slightly sub-Keplerian at most radii. 

The radial run of flow quantities in the inner regions 
can be fit by power laws, as done by I A and SPB. Our 
best fit power laws for the fiducial model (Run A) over 
2.7GM/c 2 < r < 20GM/c 2 are p cx r" 6 , c s cx r" ' 5 , 
\v r \ cx r~ 2 , and v<f, cx r~°' 8 . Between 2.7GM/c 2 < r < 
6GM/c 2 , is best fit by cx r~ 0,9 , which is nearly, but 
not exactly, consistent with conservation of fluid specific 
angular momentum (v<f, cx r -1 ). Angular momentum is 
not exactly conserved at r < 6GM/c 2 because of viscous 
torques. 

The careful reader may notice that the power law slopes 
quoted in the last paragraph are not consistent with a con- 
stant mass accretion rate. This is because the power laws 
are derived from averages over \6 — tt/2\ < ir/6, following 
IA. If one averages over all 8 and time, then the result- 
ing profiles are consistent with constant mass, energy, and 
angular momentum accretion rates, as they must be for a 
flow that is steady when averaged over large times. 

4.2. Dependence on Inner Boundary Location and 
Gravitational Potential 

have chosen to use IA's range in 8. 
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Having established that the differences between SPB 
and IA's models are due to model choices rather than nu- 
merics, we were also interested in studying whether any 
features of global viscous accretion models are strongly 
dependent on numerical parameters. The first parameter 
we considered was the location of the inner boundary. 

In models that use a Newtonian gravitational potential 
(such as IA and SPB) the location of the inner boundary 
is not an interesting parameter in the sense that there is 
no physical lengthscale that one can compare Ti n to: it is 
simply a scaling parameter. In models that use a pseudo- 
Newtonian potential, however, there is a feature (a "pit") 
in the potential on a lengthscale GM/c 2 . Starting with 
our fiducial model, then, what is the effect of shifting rj„? 

Our fiducial run has rj„ = 2.7GM / 'c 2 . This may be 
compared with Run B, which has r in = 6GM/c 2 . Figure 
A4 compares the accretion rates in the two runs. Evi- 
dently there are two changes in the solution. First, the 
time-averaged accretion rates differ by a large factor. The 
mean mass accretion rate is factor of 3 lower in Run B 
than Run A. The reduced eigenvalues e and I also differ 
by about 50% (see Table 2). Second, the time variation of 
the accretion rates differs, with Run B showing far more 
short-timescale variations. The short-timescale variations 
are due to the interaction of unstable convective modes 
with the boundary conditions. Inspection of the runs re- 
veals an enhancement of convection and turbulence near 
ri n in Run B. 

The differences between Run B and Run A are caused by 
the boundary location. Gradual variation of r in (in models 
not discussed in detail here) reveals that if the flow on the 
inner boundary is everywhere and always supersonic, then 
the solution is similar to Run A. If the flow is subsonic, 
then the solution exhibits artifacts like those seen in Run 
B. 

Evidently forcing the flow to be supersonic on the inner 
boundary causally disconnects the flow from the bound- 
ary. 4 This eliminates nonphysical reflection of linear and 
nonlinear waves from the boundary and renders the pre- 
cise implementation of the numerical boundary conditions 
irrelevant. 

We do not want the reader to think that this problem 
arises because we happened to choose the wrong numerical 
implementation of the boundary conditions. Our imple- 
mentation is the standard ZEUS outflow boundary con- 
dition, and it is widely used in astrophysical problems. 
While it may be possible to implement more transpar- 
ent boundary conditions in the context of other numerical 
schemes, a survey of the numerical literature shows that in 
multiple dimensions this is an area of active research (Roc 
1989; Kami 1991; Dedner et al. 2001; Bruneau & Creuse 
2001), and that no general solution to the problem has 
been found. 

Furthermore, a simple example shows that no local ex- 
trapolation scheme can work for all accretion problems. 
Consider a numerical model of a steady spherical inflow 
(Bondi flow) in a gravitational potential \l/(r). Let us sup- 
pose that we are primarily interested in accurately measur- 
ing M. We know from the analytic solution of the problem 
that M depends on the shape of the potential everywhere 

4 Although the viscous fluid equations of motion are not hyperbolic, 
the rest of the flow, the coupling is exponentially weak. 



outside the sonic point. If we place the inner boundary 
ri n outside the sonic point and use a local extrapolation 
scheme, we won't always get the correct answer because 
the local extrapolation doesn't have any information about 
the shape of the potential between ri n and the sonic point. 
Put differently, one can't determine a global solution from 
local extrapolation at the boundary. The key point is that, 
while aspects of the solution may be accurate, M (and L 
and E) are sensitive to the boundary conditions. 

It is worth noting that the outer boundary is always 
in causal contact with the flow, but does not cause the 
same type of artifacts as the inner boundary. Experiments 
show that the flow is qualitatively insensitive to the loca- 
tion and implementation of the outer boundary condition. 
The time averaged M, however, is sensitive to both r out 
and ri n j/r out . The time averaged I and e scale out this 
mass dependence and so are qualitatively and quantita- 
tively insensitive to both r ou t and ri n j/r out . 

It is also worth noting the effects of changing the grav- 
itational potential. Run C (identical to IA99 Model 5) is 
identical to Run B except that the potential is now Newto- 
nian. It is qualitatively similar to Run B, but M is now a 
factor of 5 lower than Run A. Run C also has the property 
that I oscillates about 0.0. This is a problem if the focus 
of the simulation is measuring M or L. 

To summarize: the location of the inner radial bound- 
ary can determine the character of the flow. If the 
flow is everywhere and always supersonic (or super-fast- 
magnetosonic in MHD) on the inner boundary then 
boundary-related corruption of the flow is impossible. 
Since it is computationally expensive to place the inner 
boundary very deep in the potential (for our model, the 
time step dt <~ (?*j„ — 2GM/c 2 )), the optimal location for 
the inner boundary is just inside the radius where the ra- 
dial Mach number always exceeds 1. 

The results of IA and SPB do not focus on the time- 
dependence of the accretion values, so much of their dis- 
cussion is unaffected by their treatment of the inner radial 
boundary. As discussed below, there are small changes 
related to the appearance of outflows. 

4.3. Comparison of Torus and Injection Models 

The torus and injection methods represent sharply dif- 
ferent approaches to studying accretion flows. The equilib- 
rium torus presents a physically well-posed problem, but 
the accretion flow is transient: no steady state can be 
achieved. The injection method reaches a quasi-steady 
state, but much of the computational domain is wasted on 
evolving the injection region, which has no astrophysical 
analog: it is nonphysical. It is natural to ask whether these 
two widely used schemes for supplying mass can be made 
comparable or used to measure any of the same quantities. 

We selected two runs, E (torus) and F (injection), that 
had similar mass distributions in an evolved state. The 
torus run was studied at a time when M was close to its 
maximum. Run F's M is a factor of 10 larger than Run 
E's. This difference might have been anticipated from the 
sensitivity of the injection run to r out and r in j/r out : runs 
in which mass is concentrated closer to the outer bound- 
ary tend to have lower accretion rates because more of the 

the flow in the supersonic region is in principle in causal contact with 
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mass escapes through the outer boundary. The time aver- 
aged mass accretion rate is therefore strongly dependent 
on the method of mass supply. 

The energy and angular momentum accretion rates per 
unit mass are, however, insensitive to the experimental de- 
sign. We find that e and I differ by less than 3% in Runs 
E and F (see Table Al and Table A2). These quantities 
are apparently set by conditions near the inner boundary 
(the ISCO for the Pseudo-Newtonian potential), and can 
be measured in either type of experiment. 

4.4. Other Parameters 

We have varied ri n j and r out /r in j, and as reported 
above, these strongly affect the time-averaged value of M. 
The sense of the effect is that a simulation with a larger 
rout/rinj loses less matter through the outer boundary, 
and this results in more matter streaming back into the 
black hole (by up to a factor of 10). The qualitative nature 
of the flow, however, is roughly independent of r out /ri n j 
in that, e.g., the temporal power spectrum of M is sim- 
ilar. The qualitative nature of the flow is dependent on 
rtnj. If one fixes r out /ri n j and all other parameters, the 
range in a where unipolar outflows are observed tends to 
become smaller and disappears altogether for r in j as small 
as AQGM/c 2 . 

The dependence of accretion models similar to ours on 
a has already been investigated by 1A. They find that the 
flow changes from turbulent to laminar as a is increased 
and the higher viscosity damps modes of increasing length- 
scale. IA find that a < 0.03 the flow is turbulent, and for 
a > 0.3 the flow is laminar. For 0.03 < a < 0.3 the flow ex- 
hibits a "unipolar" outflow. Our results are in agreement 
with IA. However, our models with a pseudo-Newtonian 
potential and super-sonic flow at the inner radial boundary 
show a slight shift in the values of a that exhibit unipolar 
outflows. 

We did find a critical value of a « 0.5 above which 
no supersonic flow at r,-„ could be achieved due to vis- 
cous heating, at least for 7 = 3/2 and 7 = 5/3 and for 
r in > 2.1GM/c 2 . Smaller values of r in were not computa- 
tionally practical. This high a is typically associated with 
a bipolar outflow, as seen by IA. Even in this case, how- 
ever, a choice of rj„ = 2.1GM/c 2 instead of rj„ = 6GM/c 2 
leads to a qualitatively different profile for the flow. The 
flow with smaller ri n — 2.1GM/c 2 has a bipolar outflow 
starting at larger radius (lQGM/c 2 ) rather than immedi- 
ately on the boundary as with rj„ = 6GM/c 2 , and the 
mass accretion rate increases by a factor of 3. 

Finally, we studied the dependence of the results on nu- 
merical resolution. We find that N r x Ng = 108 x 50 is 
sufficient at a = 0.1 to resolve the shortest wavelength 
convective mode. Also, we chose our value of r to agree 
with SPB99's torus models. We experimented with vary- 
ing r and find that, all things being equal, smaller r gives 
more laminar flow. 

5. SUMMARY 

Work in this field will shortly focus on global MHD mod- 
els in pseudo-Newtonian potentials and in full general rel- 
ativity. In our view it is useful to understand the solution 
space for physically and numerically simpler viscous mod- 
els before turning to MHD. It is even possible, as Stone & 



Pringlc (2001) have claimed, that viscous hydrodynamics 
provides a crude approximation to the MHD results. In 
any event, this investigation provides a preview of some of 
the experimental issues that will play a role in most future 
global numerical investigations of accretion flows. 

This investigation was initially motivated by a desire to 
understand whether the differences between earlier global 
viscous hydrodynamics simulations performed by IA and 
SPB were caused by differences in experimental design or 
numerical method. IA and SPB reported different degrees 
of convective turbulence in their models and found differ- 
ent radial scalings for vertically averaged quantities such 
as temperature and density. Using a single code, we were 
able to reproduce both sets of results. We conclude that 
the differences are due to experimental design. 

We also found, while reproducing IA and SPB's results, 
that some aspects of our solutions were sensitive to the nu- 
merical treatment of the region close to the inner bound- 
ary in models that use a pseudo-Newtonian potential. In 
particular, / and e, the specific angular momentum and en- 
ergy of accreted material, are strongly dependent on r in , 
the location of the inner boundary. When the flow is su- 
personic at ri n the location of the boundary does not affect 
I and e. But when the flow is subsonic at n„ the flow in- 
teracts strongly with the numerical boundary condition. 
This produces spurious outflow events and makes I and 
e dependent on rj n . Evidently for accurate measurement 
of these quantities it is necessary to isolate the numerical 
boundary condition behind a sonic transition that is lo- 
cated within the computational domain; one must place 
the inner boundary condition inside a "sound horizon" . 

We are not saying that all models that lack a sonic 
transition in the computational domain are fatally flawed. 
Whether the treatment of the inner boundary condition 
is problematic or not depends on what is being measured. 
For the nonlinear eigenvalues L, E, and M that we have fo- 
cused on here, however, the treatment of the inner bound- 
ary condition is crucial. Furthermore, the only guarantee 
that the inner boundary condition is not governing the so- 
lution is to isolate it behind a sonic transition; this is the 
only completely safe choice. 

The location of the inner boundary may prove even more 
crucial in MHD models. Much of the character of the flow 
is determined in the turbulent, energetically important re- 
gion of the flow just outside the fast magnctosonic transi- 
tion, just as the region immediately outside the sonic tran- 
sition determines the nature of the viscous flows described 
in this paper. We have performed some preliminary nu- 
merical tests and find that, as in the viscous flow, I and 
e for an MHD flow are sensitive to the treatment of the 
inner boundary. For various reasons it may prove difficult 
to achieve a fast magnetosonic transition in the computa- 
tional domain; accurate treatment of the inner boundary 
condition may require fully (general) rclativistic MHD. 

We have also compared two commonly used experimen- 
tal designs for black hole accretion flow studies: models 
that begin with an equilibrium torus, and models that 
continuously inject fluid onto the grid. The choice between 
these models is to some degree a matter of taste. We find 
the equilibrium torus slightly easier to initialize and an- 
alyze. Remarkably, the two different approaches produce 
indistinguishable measurements for I and e, the specific 
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angular momentum and energy of the accreted material. 

A parallel, viscous, axisymmetric hydrodynamics code 
based on that used in this paper can be found at 
http://kcrr.physics.uiuc.edu. 

This work was supported in part by a GE fellowship 
and NASA GSRP Fellowship Grant NGT5-50343 to JCM, 



Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214 

Balbus, S. A. & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 

Blandford, R. & McKee, C. 1982, ApJ, 255, 419 

Blandford, R. D. & Payne, D. G., MNRAS, 199, 883 

Blondin, J., Borkowski, K., & Reynolds, S. 2001, ApJ, 557, 782 

Bondi, H, 1952, MNRAS, 112, 195 

Bruneau C. H. & Creuse E. 2001, Int. J. Numer. Meth. Fluids, 36, 
807 

Dedner A. ct al. 2001, J. Comp. Phys., 171, 448 
Doeleman, S. S. et al. 2001, ApJ, 121, 2610 
Gcbhardt, K. et al.2001, ApJ, 543, 5 

Hawley, J. F., Wilson, J., & Smarr, L. 1984, ApJS, 55, 211 

Hawley, J. F. 1991, ApJ, 381, 496 

Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 223 

Hawley, J. F., Gammie, C. F.,& Balbus, S. A. 1995, ApJ, 440, 742 

Hawley, J. F. 2000, ApJ, 528, 462 

Hawley, J. F. 2001, ApJ, 554, 534 

Hawley, J. F. & Krolik, J. 2001, ApJ, 548, 348. 

Hawley, J. F., Balbus, S.A., & Stone, J. 2001, ApJ, 554, 49 

Igumenshchev, I. V. & Abramowicz, M. A. 1999, MNRAS, 303, 309 

Igumenshchev, I. V. & Abramowicz, M. A. 2000, ApJS, 130, 463 



an NCSA Faculty Fellowship for CFG, the UIUC Research 
Board, and NSF Grant AST-0093091. Computations were 
done in part under National Computational Science Al- 
liance grants AST010012N and AST010009N using the 
Origin 2000 and posic Linux cluster at NCSA. We thank 
Stu Shapiro and the referee for comments that substan- 
tially improved this paper. 



(Al) 



(A2) 



(A3) 



(A4) 



(A5) 



(A6) 



Igumenshchev, I. V., Abramowicz, M. A., & Narayan, R. 2000, ApJ, 
537, L27 

Junor, W. & Biretta, J. & Livio, M. 1999, Nature, 401, 891 
Kami, S. 1991, Int. J. Numer. Meth. Fluids, 13, 201 
Koidc, S., Shibata, K., & Kudoh, T. 1999, ApJ, 522, 727 
Koidc, S., Meier, D. L., Shibata, K., & Kudoh, T. 2000, ApJ, 536, 
668 

Lo, K. Y., Shen, Z., Zhao, J., & Ho, P. T. P. 1998, ApJ, 508, 61 
Lynden-Bell, D. & Pringlc, J. E. 1974, MNRAS, 168, 603 
Meier, D., Koidc, S., & Uchida, Y. 2001, Science, 291, 84 
Narayan, R. & Popham, R. 1993, Nature, 362, 820 
Paczynski, B. & Wiita, P. J. 1980, A&A, 88, 23 
Papaloizou, J. & Pringle, J. 1984, MNRAS, 208, 721 
Rees, M. J. 2001, preprint (astro-ph/0101266) 
Roe, P. L. 1989, Comp. Fluids, 17, 221 
Salpeter, E. E. 1964, ApJ, 140, 796 
Shakura, N. I. & Sunyacv, R. A. 1973, A&A, 24, 337 
Stone, J. M., Pringle, J. E., & Begelman, M. C. 1999, MNRAS, 310, 
1002 

Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 753 
Stone, J. M. & Pringle, J. E. 2001, MNRAS, 322, 461 
Zel'dovich Ya. B. 1964, Sov. Phys.-Dokl., 9, 195 



APPENDIX 
RATE OF STRAIN TENSOR 

The rate of strain tensor e is a symmetric tensor that in spherical polar coordinates has 

rr dv r 1 
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Fig. Al. — Time-averaged spatial structure of fiducial run (Run A; a = 0.1, r ln = 2.7GM/c 2 , and r ou t = 600GM/c 2 ). Shown arc the 
density (upper left), Bernoulli parameter (Be = (l/2)v + c 2 /(7 — 1) + W) (upper right; dotted line is a negative contour), scaled mass flux 
r 2 sinS(pv) (lower left), and scaled angular momentum flux r 3 sin 2 9(pvv t f > + II ■ <j>) (lower right). The flow is not symmetric about the equator 
because the flow exhibits long timcscalc antisymmetric variations. Convective bubbles form at the interface between positive and negative 
Bernoulli parameter (i.e. unbound and bound matter). 
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Fig. A2. — The evolution of M/M in j, e = E/(Mc 2 ), and / = L c/(GMM) in the fiducial run (Run A). The dotted line indicates the 
thin disk value. The run has clearly entered a quasi-steady state. The evolution is relatively smooth with a small variation on a timescale 
r Ri 4 X 10 4 . This is the timescale for convective bubble formation (the low point in mass accretion rate is when bubble forms). For this 
model the bubble forms at alternate poles. A full cycle requires of order one rotation period at the injection radius. 
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Fig. A3. — The radial run of 9 and time averaged quantities from the fiducial run (Run A). Shown are the density (upper left), squared 
sound speed (upper right), radial velocity (lower left), specific angular momentum (lower right; solid line), and circular orbit specific angular 
momentum (lower right; dashed line). Crudely speaking, the inner flow is consistent with a radial power law. The best fits to a power law 
are: p cx r -0 ' 6 , c s oc r -0 ' 5 , \v r \ oc r~ 2 , and vj, oc r~ ' 8 . The plots are averaged over 8 = iz/2 ± n/6. 
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Fig. A4. — The effect of moving the inner boundary on the accretion rates of mass, angular momentum, and energy (Run A vs. Run 
B). The top panel shows M/Mi„j, the middle panel E/(Mc 2 ), and the bottom panel L c/(GMM). The solid curve is Run A, which has 
ri n = 2.7GM / c 2 . The dashed curve is Run B, which has rt n = GGM/c 2 . Evidently Run B has a different variability structure and different 
time averaged values for the accretion rates. The relatively rapid and high-amplitude variations in Run B are due to nonphysical interactions 
with the inner radial boundary. Only by ensuring a supersonic flow (as in Run A) can one avoid these nonphysical effects. 
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Table Al 



Parameter List 



Run 


JV r 


Ng 


Vise. 


Potential 


Rin/r g 


Rout/r g 


Rinj /fg 


7 


a 


t f (c 3 /GM) 


A 


108 


50 


IA 


PN 


1.35 


300 


248 


3/2 


0.1 


7.3 x 10 5 


B 


80 


50 


IA 


PN 


3 


300 


248 


3/2 


0.1 


4.4 x 10 s 


C 


80 


50 


IA 


Newt. 


3 


300 


248 


3/2 


0.1 


7.3 x 10 s 


D 


64 


40 


MG 


PN 


1.2 


76 


62 


3/2 


1.0 


3.3 x 10 3 


E 


128 


80 


MG 


PN 


1.4 


21 


17 


5/3 


0.01 


3.0 x 10 4 


F 


128 


80 


MG 


PN 


1.4 


81 


21 


5/3 


0.01 


6.8 x 10 4 



Note. — IA and MG are viscosity prescription described in equations 7-9. PN is the pseudo-Newtonian 
potential of Paczynski & Wiita (1980). Run B uses Run C as initial conditions. r g — 2GM/c 2 . Ri n j for Run F 
is the position of the torus density peak po- 
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Table A2 



Results List 



Run 


Steady State Time (GM/c 3 ) 


Max. Mach at R in 


M 


-(E/(Mc 2 )) x 10 -2 


L c/(GMM) 


A 


4.3 x 10 5 


-1.4 


5.96 x 10~ 2 


2.06 


1.75 


B 


2.4 x 10 s 


+0.0 


1.95 x 10~ 2 


3.72 


1.29 


C 


2.4 x 10 s 


+0.0 


9.46 x 10~ 3 


6.77 


.0746 


D 


> 3.3 x 10 3 


-0.4 


> 3.59 x 10~ 2 


4.64 


-0.167 


E 


5.5 x 10 3 


-3.0 


3.48 x 10~ 2 


3.01 


3.41 


F 


2.0 x 10 4 


-3.4 


5.03 x 10 _1 


3.11 


3.35 



Note. — Runs A-E are injection runs with mass accretion rate units in Mi n j and Run F is a torus run with mass accretion 
rate unit in po(GM) 2 /c 3 . Run C's angular momentum fluctuations are 10 times the average value shown. 



